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A quantum Monte Carlo method combining update of the loop algorithm with the global flip 
of the world line is proposed as an efficient method for studying the magnetization process in an 
external field, which has been difficult because of inefficiency of the update of the total magneti- 
zation. The method is demonstrated in the one dimensional antiferromagnetic Heisenberg model 
and the trimer model. We attempted various other Monte Carlo algorithms to study systems in 
the external field and compared their efficiency. 
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§1. Introduction 

The quantum Monte Carlo method has become one of 
the most powerful methods to study strongly correlated 
quantum systems. Since this method was proposed, vari- 
ous improvements have been introduced. 1 ' Several prob- 
lems have been pointed out concerning the efficiency of 
the method in the original style which is called the world 
line quantum Monte Carlo method (WLQMC), such as 
long autocorrelation in Monte Carlo update, the nuisance 
of the extrapolation of the Trotter number, inefficient 
sampling in the study of the magnetization process in an 
external field 2 ' 3 ' and the negative sign problem. 4 ' 

Recently what is called the loop algorithm has been 
introduced 5 ' to overcome the problem of long autocor- 
relation in Monte Carlo update. Furthermore instead of 
discrete time with finite Trotter number, an algorithm 
using continuous time has been introduced 6 ' and the 
nuisance of the extrapolation of the Trotter number has 
been overcome. These improvements allow us to study 
systems at very low temperature. 

So far some successful approaches to the magnetiza- 
tion process have been reported. 7 ' 8 ' In this paper we 
propose a direct application of the loop algorithm with 
continuous time method (LCQMC) to the study of the 
magnetization process. There are several realizations of 
LCQMC in the external field. We attempted various 
methods and compared their efficiency. We found that 
a hybridization of the standard LCQMC and the global 
flip in WLQMC is particularly efficient for all values of 
the field. We demonstrate this method for the one di- 
mensional antiferromagnetic Heisenberg model (1DAFH) 
and also for the trimer model. 9, 10 ' 

§2. Model and Method 

In the quantum Monte Carlo method using the Suzuki- 
Trotter decomposition: 

. 



we express a d-dimensional quantum state in a (d + 1)- 
dimensional classical configuration. The new axis is 
called the Trotter direction which comes from the decom- 
position. 1 ' The thermal average of quantities is obtained 
by sampling in the classical configuration. 

In WLQMC, in order to change the total magnetiza- 
tion, what is called the global flip is carried out where a 
straight world line is flipped simultaneously. It is, how- 
ever, very rare that the global flip is accepted at low 
temperature in systems with the strong quantum fluctu- 
ation, because the straight world line exists with a very 
low probability. 3 ' On the other hand, in LCQMC, we 
generate various shapes of the loop which has nonzero 
magnetization. The flip of such loops causes change of 
the total magnetization. Thus it would be expected that 
the problem of inefficient sampling in Monte Carlo study 
of the magnetization process in an external field can be 
solved in LCQMC. 

In order to investigate the feasibility of the method, we 
study the magnetization process of the one dimensional 
antiferromagnetic Heisenberg model: 



L 

Si+i - h sf 



(2) 



where Sf = \af (a a is the Pauli matrix), a = x,y, 
and z. Here a periodic boundary condition is adopted. 
According to the standard loop algorithm, 5 ' a graph 
is assigned to each plaquette where four spins inter- 
act. The spin configurations of the plaquette and the 
types of the graph are depicted in Figs. 1(a) and 1(b), 
respectively. In Table I the weights {w g (C)} of the 
graphs {g=Gl,G2,G3,G4} in the spin configurations 
{C=C1,C2,C3,C4} for H = are listed, where 
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there 



At = 



I 

N' 



(4) 



with the Trotter number N. W(C) is the Boltzmann 
factor of the configuration C and 



W(C) = ^> g (C). 



(5) 



Graphs are allocated with the weights {w g (C)} and they 
form loops. 

In the continuous time algorithm, 6 ) 



N 



oo, 



that is ^~ 
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0, 



(6) 



a horizontal cut which corresponds to G2 in Fig. 1(b) is 
allocated with a probability density J/2 where the spins 
in both sides are antiparallel which corresponds to C2 in 
Fig. 1(a). 

The effect of the external field H is taken into account 
through the global weight. 5 ) That is, we flip a loop with 
a magnetization m with a probability 



p(m) 



exp(— f3Hm) 



exp(pHm) + cxp(—/3Hm) ' 
where m is defined by 
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Here 



fioop denotes an integration along a loop and m 
takes values 0, ±|, ±1, • • •. 

The Monte Carlo results for L = 12 are compared with 
the exact ones obtained by the diagonalization method. 
It should be noted that the efficiency of Monte Carlo 
sampling for the magnetization process becomes better 
if the length of the chain increases, because the efficiency 
depends on how often processes changing the total mag- 
netization are accepted. Generally if we can produce 
a good result for the magnetization process in a short 
chain, we expect successful results in longer chains. Thus 
in this paper we investigate the efficiency in a short chain 
L = 12. The magnetization processes were obtained with 
1,000,000 Monte Carlo steps (MCS). The simulation of 
10 6 MCS is divided into 10 bins (each bin has 10 5 MCS). 
The errorbar is estimated from the standard deviation 
of the distribution of the data of the bins. It should be 
noted that errorbar in this definition becomes very small 
if the configuration freezes in some local stable configu- 
ration during Monte Carlo simulation, where all the bins 
give almost the same value. The true errorbar should 
be obtained from the distribution of the data of several 
independent simulations. In this paper, however, the er- 
rorbar denotes the one in the former definition as far as 
no particular attention is payed. 

In Fig. 2(a) we show the magnetization process at a 
high temperature (T/J = 1.0), where the circle denotes 
the Monte Carlo result and the solid line denotes the ex- 
act value. There we find full agreement. Furthermore it 
shows excellent convergence of the data with small error. 
Thus the method produces good results as was expected 
at high temperature. In Fig. 2(b) we show the magneti- 



zation process at a low temperature (T/J = 0.1). There 
the present method reproduces the exact data at low 
fields up to H ~ J. At higher fields, however, the method 
fails to produce the correct magnetization, where a kind 
of freezing of the configuration seems to occur. This ob- 
servation can be interpreted as follows. The loops to be 
flipped are made by assigning graphs with the weight 
in Table I which is suitable for H = 0. Thus update 
with the flip of such loops is considered to be an im- 
portant sampling for H = and efficient to study low 
field properties. On the other hand, this sampling is no 
more efficient at high field, although it would provide the 
correct result with infinite MCS in principle. The config- 
uration remains in a fixed magnetization and the method 
can not produce the correct result. In order to overcome 
this difficulty, we have to introduce other recipes. In 
the following, we examine two methods to overcome this 
difficulty. 

2.1 LCQMC using weights for H ^ 

One plausible way to overcome the difficulty is to use 
graphs for H ^ 0. The weight w e (C) consists of two 
parts: 



Wg (C) = % A g (C), 



(9) 



where A g (C) = 1 or 0. The weight for H = shown in 
Table I corresponds to a set of {A g (C)} shown in Table 
II and 
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in the first order of At. For H ^ the Boltzmann 
factors of CI and CI' are different. In this case there is 
no solution for the graph weights {v g } within {A g (C)} 
given in Table II. Thus we have to choose another set 
of {Ag(C)} and look for a positive solution {v g } for the 
eqs. (9) and (5). For example if we choose a set shown 
in Table III, we find a solution 

At At 

wgi = 1 - — J - —H, 

At 

VG2 = —J, 
«G3 = — H, 
«G4 = — H. 

Because the effect of the field is taken into account in the 
graph weight, we flip each loop with a probability 1/2. 
However, the new weight violates the detailed balance. 
In order to recover the detailed balance, some types of 
flip must be prohibited. In the above case we do not flip 
loops which have either the cross cut G3 or the freezing 
cut G4 in the configuration CI. 

In Fig. 3(a) we show the magnetization process at a 
high temperature (T/J = 1.0), where the circle denotes 
the Monte Carlo result and the solid line denotes the ex- 
act value. There we find full agreement as well as the 
standard LCQMC. In Fig. 3(b) we show the magneti- 
zation process at a low temperature (T/J = 0.1). The 
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data are consistent with the exact value at low fields. 
At high fields, however, it seems that a kind of freezing 
of the configuration occurs as was observed in the stan- 
dard LCQMC. Change of the total magnetization hardly 
occurs at high field and low temperature, because the 
number of loops with prohibited cuts increases. 

There are several other choices of {A g (C)} where we 
can find a positive solution {v g }. We also examined such 
choices but we found the similar difficulty in all the cases. 

2.2 Hybridization of LCQMC and the global flip in 
WLQMC 

In this circumstance we introduce the global flip of the 
world line, remembering WLQMC can produce M(H) at 
high field rather well. We perform the following pro- 
cedures. First we perform one step of the standard 
LCQMC and then we look for a straight world line and 
flip the whole line with a probability 



Pgf = 



w(-S z ) 



(12) 



where S z is the ^-component of the spin in the straight 
line and w(S z ) is defined as 

\og(w(S z )) = V / /J Iog(W J -(r))dr + ) 9H5* ) (13) 
i J ° 

where the summation over j denotes all bonds connected 
to the concerning site, i.e., here the left and right near- 
est neighbors. When the spin in the neighboring line is 
parallel or antiparallel, the weight Wj(r) is given by 



log(W,(r)) = --, 



or 



log(W j (r)) = 



4' 



(14) 



(15) 



respectively. 

In Fig. 4(a) we show the magnetization process at a 
high temperature (T/J = 1.0), where the circle denotes 
the Monte Carlo result and the solid line denotes the ex- 
act value. There we find full agreement as well as the 
previous methods. In Fig. 4(b) we show the magnetiza- 
tion process at a low temperature {T/J = 0.1). Here we 
find a good result even at the low temperature. There 
the uneven magnetization process due to the finiteness of 
the chain is also reproduced correctly with small error. 
We find this hybrid method practically useful, although 
the method has also the difficulty at much low tempera- 
ture. 

Here it should be noted that if we start from a fully 
magnetized state, a straight world line appears as a loop 
in LCQMC. Thus one may consider that M(H) can be 
correctly produced only by LCQMC if we start from a 
high field with the fully magnetized state and then re- 
duce the field gradually. So far we started simulation 
for the case of H = and then continued simulation 
increasing the field gradually without reset of the spin 
configuration. In Fig. 5 we show results obtained by 
simulation starting from the high initial field for both the 
standard LCQMC and the hybrid method. The reason 



of this failure in the standard LCQMC is the following. 
In LCQMC, the probability of the flip of a straight world 
line surrounded by parallel straight lines is 

1 



PLC = c /»r + r (16) 

while the probability of the flip of the line in the global 
flip is 

Thus we find that the escape rate from the metastable 
configuration, i.e., a flip of the straight world line sur- 
rounded by parallel ones, is much larger in the global 
flip. Thus the hybrid method gives better result. 

§3. Magnetization Process of the TVimer Sys- 
tem 

The compound 3CuCi2-2dioxane has been studied as a 
trimer system. 9 ) The magnetization process of the trimcr 
system has attracted interest. 10 - 1 The difficulty of Monte 
Carlo simulation of WLQMC was pointed out. 2 ' That 
is, if we use a large Trotter number, change of the total 
magnetization hardly occurs as we mentioned above. In 
such situation, if we try to estimate the extrapolated 
value from rather small values of the Trotter number, it 
is quite possible to conclude an apparent extrapolated 
value which is very different from the true one. 

Here we apply the hybrid method to this model. The 
Hamiltonian of the model is given by 



H = 



L-2 

E [ 

-l,step3 



+JAFS i+2 ■ Si + s] - H S Z . 

First we study the magnetization process for 
7 = Jf/Jaf = 5, 



(18) 



(19) 



which is close to the ratio for the compound 3CuCi2 • 
2dioxane. 2 ' 9 -' In Fig. 6 we show M(H) for L = 12 at a low 
temperature (T/Jaf = 0.1), where the circle denotes the 
Monte Carlo result and the solid line denotes the exact 
value. There we find full agreement. We also found that 
we can obtain good results even by the standard LCQMC 
in the trimcr system. The reason is that the present 
system contains ferromagnetic bonds where we allocate 
the cross cut G3. This cut causes the loop to prolong 
along the Trotter axis and thus the number of loops with 
nonzero magnetization increases. The flip of such loops 
contributes to the equilibration of the magnetization. 

Next we study temperature dependence of M(H) for 
7 = 4.4, which has been estimated for the compound 
3CuCl 2 • 2dioxane by Hida (J F = 136 K, J AF = 30.6 
K). 10 ) In Fig. 7 we show the Monte Carlo data of M(H) 
for L = 60, which is considered to be long enough to 
represent the magnetization process in the thermody- 
namic limit. We find a similar temperature dependence 
of M(H) to the experimental data. However, M(H) at 
low temperatures (T/Jaf = 0.049 and 0.14) in the ex- 
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periment show steeper gradient near the saturated field 
than those in Fig. 7. This disagreement suggests that 
the Hamiltonian (18) does not represent the experimen- 
tal situation completely, although it gives a good approx- 
imation of the system. Attempt at more precise tuning 
of the parameters to fit the data would provide further 
information of the material, which is an interesting prob- 
lem in the future. 

When we decrease the ratio 7, the system shows a 
1/3 plateau in the magnetization process. 10 ) In Fig. 8 
we show temperature dependence of M(H) for the ratio 
7 = 1 for L = 60. There we find the plateau at very low 
temperature (T/Jaf = 0.05). The step like structure is 
smoothed out at rather low temperature (T/Jaf = 0.1) 
and a smooth curve is obtained at a high temperature 
(T/Jaf = 1.0). 

At low temperature and high field the autocorrelation 
of the simulation is large even in the present improved 
method and the errorbar defined previously does not rep- 
resent the correct variation of the data. For such cases 
we performed five independent runs and estimated the 
standard deviation of data over the runs, which indicates 
how much data scatter. In Fig. 8 the errorbars for data 
shown by painted symbols denote the standard deviation 
of the data in five independent runs of 10 6 MCS, while 
the errorbars for data shown by open symbols denote 
the errorbars in a run of 10 6 MCS. At high tempera- 
tures (T/Jaf = 0.1 and 1.0), the errorbars estimated 
from independent runs are small, where the single run is 
considered to be enough. 

§4. Summary and Discussion 

We attempted several methods to find efficient Monte 
Carlo algorithm for the magnetization process and found 
a hybridization of the standard LCQMC and the global 
flip in WLQMC to be efficient and practically useful. 
We applied the method to the one dimensional antifer- 
romagnetic Heisenberg model which is a lattice consist- 
ing of only antiferromagnetic bonds, which is regarded as 
the hardest model to study the magnetization process by 
quantum Monte Carlo methods. Using the method we 
also obtained the temperature dependence of the magne- 
tization process of the trimer system successfully, which 
has been difficult in WLQMC. 

So far various interesting quantum phases have been 
proposed in the ground state. However, usually in exper- 
iments only data at finite temperatures are available. If 
we see the properties at finite temperature at zero field, 
it is rather difficult to find characteristic properties of the 
quantum phases. Thus the data at finite field are use- 
ful to grasp the characteristics. The present method is 
powerful for such purposes. In particular the position of 
the cross of the magnetization processes is a point below 
which the magnetization decreases as the temperature 
decreases while above which it increases as the temper- 
ature decreases. Thus how the magnetization processes 
cross each other may represent a characteristic of the 
system. It would be an interesting problem to charac- 
terize properties at finite temperature for various types 
of quantum phases even if the phase transition itself is 
defined only in the ground state. 



Among the methods 7 ' 8 ) including the present method, 
a more detailed study will be required to clarify in what 
situation a particular approach is the most efficient. 
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Fig. 1. (a) Spin configurations of the plaquette and (b) types of 
the graph 



Fig. 2. Magnetization process of the 1DAFH model obtained by 
the standard LCQMC with 10 6 MCS. (a) T/J = 1.0 and (b) 
T/J = 0.1. The solid line denotes the exact values obtained by 
the diagonalization method. 



Fig. 3. Magnetization process of the 1DAFH model obtained by 
LCQMC with the Table III with 10 6 MCS. (a) T/J = 1.0 and 
(b) T/J = 0.1. The solid line denotes the exact values obtained 
by the diagonalization method. 



Fig. 4. Magnetization process of the 1DAFH model obtained by 
the hybrid method with 10 6 MCS. (a) T/J = 1.0 and (b) T/J = 
0.1. The solid line denotes the exact values obtained by the 
diagonalization method. 



Fig. 5. Magnetization process of the 1DAFH model obtained by 
the standard LCQMC (o) and the hybrid method (•) starting 
from a full magnetized state with 10 6 MCS. T/J = 0.1. The solid 
line denotes the exact values obtained by the diagonalization 
method. 



Fig. 8. Magnetization process of the trimer model obtained by 
the hybrid method with 10 6 MCS. 7 = 1 and L = 60. The open 
circle, diamond and square denote T/Jap = 0.05, 0.1 and 1.0, 
respectively. The solid symbols denote the data obtained by five 
independent runs of 10 6 MCS. 



Fig. 6. Magnetization process of the trimer model obtained by 
the hybrid method with 10 6 MCS. 7 = 5 and T/J A p = 0.1. The 
solid line denotes the exact values obtained by the diagonaliza- 
tion method. 



Fig. 7. Magnetization process of the trimer model obtained by 
the hybrid method with 10 6 MCS. 7 = 4.4 and L = 60. The sym- 
bols circle, diamond, square and triangle denote T/Jaf = 0.049, 
0.14, 0.23 and 0.29, respectively. These temperatures correspond 
to T = 1.5 K, 4.2 K, 7.0 K and 9.0 K, respectively. The field 
H = 1.0 corresponds to 21 Tesla. 
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Fig. 4 (a) 
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Table I weight for H = 
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Table II A g (C) for H = 
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